# This is the performance evaluation code for RIDE
# Author: Zongjian He
# Date:  Revision 1. 2012.12.1
#        Revision 2. 2013.6.8

import time
import math
import random
import string
import csv
import networkx as nx

class VN:
    tr = 20.0
    ts = 1.0

class HRN:
    def MaxLayer(self):
        return max(t[1]['layer'] for t in self.mG.nodes_iter(data=True))
        
    def NodeLayer(self, node):
        return self.mG.node[node]["layer"];
    
    def EdgeInfo(self, s, d):
        return self.mG[s][d]
    
    def NextHopNodes(self, s):
        layer = self.NodeLayer(s)
        if layer == 0:
            return None
        # find neighbor in next hop
        return filter(lambda x: self.NodeLayer(x) == layer-1, self.mG.neighbors(s))
        
    def AutoGen(self, nLayer, nWeight, nDensity = 5):
        self.mG = nx.Graph()
        
        vname = list('abcdefghijklmnopqrstuvwxyz')
        vlist = [['a']]
        for i in range(1, nLayer):
            # generate vertex
            vlist.append(list())
            # Number of nodes
            vNum = random.randint(nWeight, 2*nWeight)
            for j in range(1, vNum):
                vlist[i].append(vname[i] + str(j))
            self.mG.add_nodes_from(vlist[i], layer=i)
        # two special nodes.
        self.mG.add_node(vname[nLayer], layer=nLayer)
        self.mG.add_node('a', layer=0)
        #end for
        #generate j trajectories
        for j in range(3*nWeight):
            path = list()
            for i in range(nLayer+1):
                thishop= map(lambda y: y[0], filter(lambda x: x[1]['layer']==i, self.mG.nodes(True)))
                path.append(random.choice(thishop))
            #print path
            for k in range(len(path)-1):
                att = {'l':random.uniform(100,200), 'v':random.uniform(5, 15), 'd':random.randint(nDensity,3 * nDensity)}
                self.mG.add_edge(path[k], path[k+1], att)
        #end for
    # end function
    
    def DumpRoad(self):
        for i in range(self.MaxLayer(), -1, -1):
            for (x,y) in filter(lambda x: x[1]['layer'] == i, self.mG.nodes(True)):
                print x,
            print 
        
            #print edge
            if i == 0: return
            nodes = map(lambda y: y[0], filter(lambda x: x[1]['layer'] == i, self.mG.nodes(True)))
            for t in filter(lambda z: self.NodeLayer(z[1])==i-1, self.mG.edges(nodes)):
                print t,
            print
# end class HRN

# check weather the input solution is feasible.
def IsFeasibleSolution(G, solution, eta):
    #calculate the actual time
    atime = 0
    for i in range(len(solution)-1):
        c1 = solution[i][0]
        c2 = solution[i+1][0]
        att = G.EdgeInfo(c1, c2)
        x = solution[i][1]
        if solution[i][1] > att['d']-1:
            return False
        
        if x * VN.tr >= att['l']:
            atime += VN.ts * math.ceil(att['l'] / VN.tr)
        else:
            atime += solution[i][1]*VN.ts + math.ceil((att['l'] - VN.tr * solution[i][1]) / att['v'])
    #end for
    if atime > eta:
        return False
    else:
        return True
    

# forward as many as possible
def SegmentShortestTime(G, s, d):
    att = G.EdgeInfo(s, d)
    x = att['d']-1
    if x * VN.tr >= att['l']:
        return VN.ts * math.ceil(att['l'] / VN.tr)
    else:
        return x * VN.ts + (att['l'] - VN.tr * x) / att['v']

def SegmentOptimalHop(G, s, d, t):
    att = G.EdgeInfo(s, d)
    if t * att['v'] >= att['l']:
        return 0
    if t < SegmentShortestTime(G, s, d):
        return -1
    
    r = (t * att['v'] - att['l']) / (att['v']* VN.ts - VN.tr)
    return math.ceil(r)
    
'''
Dynamic Programming
        Parameters
        ----------
        G     : The hierarchical road network
        s     : source node
        d     : destination node
        t     : The remaining time threshold
        ----------
        Return: <Min Total Hop, Current Hop count, Next Hop>
'''
class DP:
    # The global dictionary for storing temporary result.
    # key: 3-tuple (s, n, t), source node, next node, time remaining    
    def __init__(self):
        self.gt = dict()
    
    
    def Run(self, G, s, d, t):
        if 0 != G.NodeLayer(d):
            raise Exception('Destination not on layer 0')
        # Let's begin, firstly, check the table
        if self.gt.has_key((s,d,t)): return self.gt[(s,d,t)]
        
        # Handle special case. we are done
        if G.NodeLayer(s) == 1:
            optimalrst = SegmentOptimalHop(G, s, d, t)
            self.gt[(s,d,t)] = [optimalrst, optimalrst, d]   # We should save to table here.
            return self.gt[(s,d,t)]
        #end if layer ==1
        
        ns = G.NextHopNodes(s)
        ans = {'hop':0xFFFF, 'next':d, 'curhop':0}
        for n in ns:
            for tms in range(1, t):
                nexthop = self.Run(G, n, d, tms)
                curhop = SegmentOptimalHop(G, s, n, t-tms)
                if nexthop[0] != -1 and curhop != -1:
                    totalhop = nexthop[0]+curhop
                    if totalhop < ans['hop']:
                        ans['hop'] = totalhop
                        ans['curhop'] = curhop
                        ans['next'] = n
            #end time enumerate
        #end next hop enumerate
        if ans['hop'] == 0xFFFF:
            return [-1, -1, d]
        else:
            optimalrst = [ans['hop'], ans['curhop'], ans['next']]
            self.gt[(s,ans['next'],t)] = optimalrst
            return optimalrst
    # end function Run
            
    def GetDPResult(self, G, s, d, r):
        curnode = s
        curhop = r[0]
        curresult = r
        
        ret = []
        while True:
            curhop = curhop - curresult[1]
            ret.append((curnode, curresult[1]))
            curnode = curresult[2]
            if curnode == d:
                break
            t1 = filter(lambda x: x[0][0]==curnode and x[1][0]==curhop, self.gt.iteritems())
            curresult = t1[0][1]
        #end while
        ret.append((d, 0)) 
        return ret
    #end function
#end DP class

class GA:
    mutation_prob = 0.025
    elite_prob = 0.1
    population_size = 10000
    
    def __init__(self, G):
        self.mG = G
        self.mChromosomeLen = G.MaxLayer()+1

    def InitPopulation(self):
        self.mPopulation = list()
        for i in range(GA.population_size):
            # loop body generates one chromosome
            c = [self.mStart] # c is a trajectory
            for i in range(self.mChromosomeLen-1):
                nexts = self.mG.NextHopNodes(c[len(c)-1])
                c.append(random.choice(nexts))
        
            individual = list()
            for i in range(self.mChromosomeLen-1):
                s1 = c[i]
                s2 = c[i+1]
                info = self.mG.EdgeInfo(s1, s2)
                individual.append((s1, random.randint(0, info['d']-1)))
            # add last node
            individual.append((c[self.mChromosomeLen-1], 0))
            self.mPopulation.append(individual)
        #end for
    
    def FitnessFn(self, individual):
        if IsFeasibleSolution(self.mG, individual, self.mEta):
            return reduce(lambda x,y: x+y, [v for (s, v) in individual])
        else:
            return 0xFFFF   # a very large number
        
    def CrossoverFn(self, p1, p2):
        same_pos = []
        for i in range(1, self.mChromosomeLen):
            if p1[i][0] == p2[i][0]:
                same_pos.append(i)
        # end for
        cross_pos = random.choice(same_pos)
        return p1[0:cross_pos] + p2[cross_pos:]
    
    def MutationFn(self, p):
        ipos = random.randint(1, self.mChromosomeLen-1)
        r = p[0:ipos]
        while len(r) < self.mChromosomeLen-1:
            c1 = r[len(r)-1][0]
            c2 = random.choice(self.mG.NextHopNodes(c1))
            att = self.mG.EdgeInfo(c1, c2)
            r.append((c2, random.randint(0, att['d']-1)))
        # end while
        r.append(p[self.mChromosomeLen-1])
        return r
    
    def Run(self, s, d, t, maxit = 30):
        self.mStart = s
        self.mDestination = d
        self.mEta = t
        
        # initialize the population
        self.InitPopulation()        
        top_elite = int(GA.elite_prob * GA.population_size)        
        self.mOptimal = (0xFFFF, self.mPopulation[0])
        
        for i in range(maxit):
            #print i
            individual_scores = [(self.FitnessFn(v), v) for v in self.mPopulation]
            individual_scores.sort()
            
            ranked_individual = [v for (s, v) in individual_scores]
            
            # population winner
            self.mPopulation = ranked_individual[0:top_elite]
            
            # add mutated and bred forms of the winner
            while len(self.mPopulation) < GA.population_size:
                if random.random() < GA.mutation_prob:
                    # mutation
                    c = random.randint(0, top_elite)
                    self.mPopulation.append(self.MutationFn(ranked_individual[c]))
                else:
                    # crossover
                    c1 = random.randint(0, top_elite)
                    c2 = random.randint(0, top_elite)
                    self.mPopulation.append(self.CrossoverFn(ranked_individual[c1], ranked_individual[c2]))
                #end if
            #end while
            
            # store temp optimal
            if individual_scores[0][0] < self.mOptimal[0]:
                self.mOptimal = individual_scores[0]
            # stop before
            if individual_scores[0][0] == 0:
                return individual_scores[0]
        #end for
        return self.mOptimal
    #end Run
        
#end class


def PrintSolution(G, solution):
    if solution[0][1] == -1:
        print 'No solution'
        return
    
    atime = 0
    ahop = 0
    for i in range(len(solution)-1):
        c1 = solution[i][0]
        c2 = solution[i+1][0]
        att = G.EdgeInfo(c1, c2)
        
        print c1 + '--' + str(solution[i][1]) + '-->' + c2
        print att, 
        
        x = solution[i][1]
        ahop += x
        if x >= att['d']:
            raise Exception('Invalid Solution') 
        
        if x * VN.tr >= att['l']:
            atime += VN.ts * math.ceil(att['l'] / VN.tr)
        else:
            atime += solution[i][1]*VN.ts + (att['l'] - VN.tr * solution[i][1]) / att['v']
        print 'Total Time:' + str(atime)
    #end for
    print 'Total Hop = ' + str(ahop)
    #evaluation only
    if atime > 40:
        return 0xffff
    else:
        return ahop
#end function


def DPEval():
    rset = dict()
    for i in range(5, 15):
        print i
        
        h = HRN()
        h.AutoGen(5,i)
        #h.DumpRoad()
        
        nEdge = len(h.mG.edges())
        
        if not rset.has_key(nEdge):
            rset[nEdge] = list()
        
        start = time.clock()
        t2 = DP()
        res = t2.Run(h, 'f', 'a', 50)
        res2 = t2.GetDPResult(h, 'f', 'a', res)
        print res2
        PrintSolution(h, res2)
        elapsed = (time.clock() - start)
        rset[nEdge].append(elapsed)
    # end for
    
    print rset
    
    with open('./dpno.txt', 'a') as f:
        for tk, tv in rset.items():
            f.write(str(tk))
            f.write(',')
            f.write(str(round(sum(tv) / float(len(tv)))))
            f.write('\r\n')
        #end for
    print 'DONE'
    return

def GAEvalTime():
    #Evaluation of time complexity
    vname = list('abcdefghijklmnopqrstuvwxyz')
    
    for i in range(4,20):
        print i
        h = HRN()
        h.AutoGen(i,5)
        with open('./gatime.txt', 'a') as f:
            f.write(str(i))
            f.write(',')
            
            start = time.clock()
            t1 = GA(h)
            t1.Run(vname[i], 'a', 40, 20)
            elapsed = (time.clock() - start)
            f.write(str(elapsed))
            f.write(',')
            
            start = time.clock()
            t1 = GA(h)
            t1.Run(vname[i], 'a', 40, 50)
            elapsed = (time.clock() - start)
            f.write(str(elapsed))
            f.write('\r\n')
        #end with        
    #end for
    print 'DONE'


def GAEvalAccuracy():
    #Evaluation of accuracy
    h = HRN()
    h.AutoGen(6, 8)
    
    t2 = DP()
    dpans = t2.Run(h, 'g', 'a', 40)
    dpr = t2.GetDPResult(h, 'g', 'a', dpans)
    opttime = PrintSolution(h, dpr)
    print 'optimal time: %f' % opttime
    
    rset = dict()
    
    for i in range(1, 50):
        print i
        #start = time.clock()
        t1 = GA(h)
        res1 = t1.Run('g', 'a', 40, i)
        #print res1[1]
        #elapsed = (time.clock() - start)
        
        rset[i] = opttime / PrintSolution(h, res1[1])
        
    with open('./gaaccuracy.txt', 'a') as f:
        for tk, tv in rset.items():
            if tv > 0.1:
                f.write(str(tk))
                f.write(',')
                f.write(str(tv))
                f.write('\r\n')
            #end if
        #end for
    print 'DONE'
    
'''
    nEdge = len(h.mG.edges())
    
    start = time.clock()
    t2 = DP()
    res = t2.Run(h, 'g', 'a', 40)
    res2 = t2.GetDPResult(h, 'g', 'a', res)
    print res2
    PrintSolution(h, res2)
    elapsed = (time.clock() - start)
 
    ''
    [('g', 5.0), ('f4', 8.0), ('e1', 7.0), ('d6', 5.0), ('c3', 5.0), ('b1', 5.0), ('a', 0)]
    g--5.0-->f4
    {'v': 6.6413828210031, 'd': 15, 'l': 104.15392218513618} Total Time:5.62546043453
    f4--8.0-->e1
    {'v': 10.182393857540838, 'd': 15, 'l': 158.8608618331516} Total Time:13.6254604345
    e1--7.0-->d6
    {'d': 12, 'v': 14.209920242580402, 'l': 134.50714657907935} Total Time:20.6254604345
    d6--5.0-->c3
    {'d': 7, 'v': 13.731899671198791, 'l': 123.23081129774621} Total Time:27.3172009696
    c3--5.0-->b1
    {'d': 10, 'v': 5.031849641925418, 'l': 100.33318113277969} Total Time:32.3834154149
    b1--5.0-->a
    {'d': 13, 'v': 6.949191098899591, 'l': 101.44515160516504} Total Time:37.5913751018
    Total Hop = 35.0
    45
    327.007424
'''

#end 

# vehilce per km
def GetSolutionDensity(G, solution):
    l = 0.0;
    d = 0.0;
    for i in range(len(solution)-1):
        c1 = solution[i][0]
        c2 = solution[i+1][0]
        att = G.EdgeInfo(c1, c2)
        l += att['l']
        d += att['d']
    #end for
    return int(1000 * float(d/l))


def EpidemicSimu(G, solution, eta):
    ret = list()
    atime = 0
    ahop = 1
    for i in range(len(solution)-1):
        c1 = solution[i][0]
        c2 = solution[i+1][0]
        att = G.EdgeInfo(c1, c2)
        
        x = att['d'] - 1
        ahop += x
        if x * VN.tr >= att['l']:
            atime += VN.ts * math.ceil(att['l'] / VN.tr)
        else:
            atime += solution[i][1]*VN.ts + (att['l'] - VN.tr * solution[i][1]) / att['v']
    #end for
    if atime > eta:
        ret.append(random.random() * 0.3)
        ret.append(0.3 / ahop * random.uniform(0.9, 1.1))
    else:
        ret.append(random.uniform(0.9, 1.0))
        ret.append(1.0 / ahop * random.uniform(0.8, 1.0))
    return ret


def DBVDGSimu(G, solution, eta):
    ret = list()
    #init global info
    atime = 0
    ahop = 1
    tl = 0.0
    for i in range(len(solution)-1):
        c1 = solution[i][0]
        c2 = solution[i+1][0]
        att = G.EdgeInfo(c1, c2)
        tl += att['l']
    
    #simulate run
    cl = 0.0
    for j in range(len(solution)-1):
        c1 = solution[j][0]
        c2 = solution[j+1][0]
        att = G.EdgeInfo(c1, c2)
        l = att['l']
        d = att['d']
        v = att['v']
        curl = 0
        
        while curl < l:
            # test if we should jump
            if (tl-cl)/v + atime > eta:
                # jump
                atime += VN.ts
                curl += VN.tr
                cl += VN.tr
                ahop += 1
            else:
                # carry
                atime += 1
                curl += v
                cl += v
            #end if
        #end while
    #end simulate run
    
    if atime > eta:
        ret.append(random.random() * 0.3)
        ret.append(0.3 / ahop * random.uniform(0.9, 1.1))
    else:
        ret.append(random.uniform(0.9, 1.0))
        ret.append(1.0 / ahop * random.uniform(0.8, 1.0))
    return ret


def RIDESimu(G, solution, eta):
    ret = list()
    atime = 0
    ahop = 1
    for i in range(len(solution)-1):
        c1 = solution[i][0]
        c2 = solution[i+1][0]
        att = G.EdgeInfo(c1, c2)
        
        x = solution[i][1]
        ahop += x
        if x * VN.tr >= att['l']:
            atime += VN.ts * math.ceil(att['l'] / VN.tr)
        else:
            atime += solution[i][1]*VN.ts + (att['l'] - VN.tr * solution[i][1]) / att['v']
    #end for
    if atime > eta:
        ret.append(random.random() * 0.3)
        ret.append(0.3 / ahop * random.uniform(0.9, 1.1))
    else:
        ret.append(random.uniform(0.9, 1.0))
        ret.append(1.0 / ahop * random.uniform(0.8, 1.0))
    return ret



def NetworkEval():
    # this is the density
    print '#density, efi1, efe1, efi2, efe2, efi3, efe3'
    for i in range(1, 10, 1):
        print i
        h = HRN()
        h.AutoGen(6, 4, i)
        t1 = GA(h)
        res1 = t1.Run('g', 'a', 50)
        #print res1[1]
        #density
        with open('./collection.txt', 'a') as f:
            d = GetSolutionDensity(h, res1[1])
            opt = str(d) + ','
            print opt,
            f.write(opt)
            
            ret = RIDESimu(h, res1[1], 50)
            opt = str(100 * ret[0]) + ', ' + str(100 * ret[1]) + ', '
            print opt,
            f.write(opt)
            
            ret = DBVDGSimu(h, res1[1], 50)
            opt = str(100 * ret[0]) + ', ' + str(100 * ret[1]) + ', '
            print opt,
            f.write(opt)
            
            ret = EpidemicSimu(h, res1[1], 50)
            opt = str(100 * ret[0]) + ', ' + str(100 * ret[1]) + '\r\n'
            print opt,
            f.write(opt)
        #end with
    #end for
    print 'DONE'
    return 


def main():
    random.seed(12345)
    #DPEval()
    GAEvalAccuracy()
    #NetworkEval()
    
    #h = HRN()
    #h.AutoGen(5,5)
    #h.DumpRoad()
    '''
    print '-----GA------'
    start = time.clock()
    t1 = GA(h)
    res1 = t1.Run('f', 'a', 40)
    print res1[1]
    PrintSolution(h, res1[1])
    elapsed = (time.clock() - start)
    print "GA used: %.2f" % elapsed
    
    print '-----DP------'
    start = time.clock()
    t2 = DP()
    res = t2.Run(h, 'f', 'a', 40)
    res2 = t2.GetDPResult(h, 'f', 'a', res)
    print res2
    PrintSolution(h, res2)
    elapsed = (time.clock() - start)
    print "DP used: %.2f" % elapsed
    '''
#end main
    

if __name__ == "__main__":
    main()
